AntibodyForests vignette: a comprehensive guide for importing 10x single-cell V(D)J sequencing data, building lineage trees, and comparing trees within and across repertoires

The VDJ_build() function imports Cell Ranger output into an R dataframe, priming the data for downstream analyses. It is optimized for Cell Ranger v7 and older versions and is suitable for B and T cell repertoires. Upon execution, the function generates a VDJ dataframe where VDJ and VJ transcripts are paired (if present) in each row, with each row representing a single cell. Post VDJ dataframe construction, a table is displayed on the console, delineating the number of filtered cells, contingent upon the remove.divergent.cells and complete.cells.only parameters.

For seamless integration with functions such as the Af_build() function (discussed next), it is strongly recommended to set the remove.divergent.cells and complete.cells.only parameters to TRUE. Cells with more than one VDJ or VJ transcript, referred to as ‘divergent cells’, typically represent doublets that may arise during the single-cell sorting, or show dual expression of kappa and lambda light chains. Such occurrences are rare, accounting for only 0.2-0.5% of peripheral blood B cells, as reported by Giachino et al. (1995). Cells harboring both a single VDJ and a VJ transcript may be included for comprehensive downstream analyses.

The output of 10x Genomics contains one output folder per run/sample. In each sample folder, the following files should be present in order to run the VDJ_build() function successfully. More information on the outputs of 10x Genomics can be found on this page.

• filtered_contig_annotations.csv
  This CSV file contains the amino acid and nucleotide framework (FWRs) and complementarity determining regions (CDRs) of the contig sequences stored in the filtered_contig.fasta file. These annotations are required to obtain the trimmed V(D)J sequences.
  
• filtered_contig.fasta
  This FASTA file contains the contig sequences from barcodes that passed the filtering steps. These contigs are annotated in the filtered_contig_annotations.csv.
  
• consensus_annotations.csv
  This CSV file contains the amino acid and nucleotide framework (FWRs) and complementarity determining regions (CDRs) of the consensus sequences stored in the consensus.fasta file.
  
• consensus.fasta
  This FASTA file contains the consensus sequence of each assembled contig, which is identical to the sequence of the most frequent exact subclonotype, in which an exact subclonotype is defined as a group of cells with identical V(D)J transcripts.
  
• concat_ref.fasta
  This FASTA file contains the concatenated V(D)J reference segments for the segments detected on each consensus sequence. Note, these reference sequences serve as an approximate reference for each consensus sequence, and need to be aligned and trimmed in order to have the same length as the consensus sequence (see trim.and.align parameter).
  
  NOTE: In older versions of Cell Ranger (older than v6), the FWR1-4 and CDR1-3 regions are not annotated. In this scenario, the all_contig_annotations.json file will be employed to obtain the ‘L-REGION+V-REGION’ and ‘J-REGION’ annotations. The raw sequences can now be trimmed using the ‘L-REGION+V-REGION’ start and ‘J-REGION’ end indices. If the all_contig_annotations.json file is not found in the sample/run directory, the columns with the trimmed sequences will remain empty and a warning message is printed.
  
  

VDJ.directory

• Description: Path to the parent directory containing the output folders (one folder for each run/sample) of Cell Ranger (note: when specified, the VDJ.sample.list parameter should not be specified).
  
• Example: "/path/to/VDJ/directory"
  
• Default: None
  
  

VDJ.sample.list

• Description: Character vector of the paths to the output folders (one folder for each run/sample) of Cell Ranger (note: when specified, the VDJ.directory parameter should not be specified).
  
• Example: c("/path/to/sample1", "/path/to/sample2")
  
• Default: None
  
  

remove.divergent.cells

• Description: If TRUE, cells with more than one VDJ transcript or more than one VJ transcript will be excluded from the VDJ dataframe.
  
• Default: FALSE
  
  

complete.cells.only

• Description: If TRUE, only cells with both a VDJ transcript and a VJ transcript are included in the VDJ dataframe.
  
• Default: FALSE
  
  

trim.germlines

• Description: If TRUE, the raw germline sequence of each clone will be aligned with the consensus sequence of that clone serving as the reference sequence, utilizing the Biostrings::pairwiseAlignment() function. The alignment conducted will be of the ‘global-local’ type, therebyy trimming the raw germline sequence.
  
• Default: FALSE
  
  

fill.germline.CDR3

• Description: If TRUE, the trimmed germline sequence of each clone will be aligned with the most frequently observed sequence for the CDR3 region of that clone serving as the reference sequence, utilizing the Biostrings::pairwiseAlignment() function. The alignment conducted will be of the ‘global-local’ type, thereby selecting the CDR3 region from the trimmed germline sequence. After the alignment, the germline CDR3 sequence is replaced by the most frequently observed CDR3 sequence, in order to obtain germline sequences that are more likely to encode producible and productive antibodies.
  
• Default: FALSE
  
  

gap.opening.cost

• Description: The cost for opening a gap when aligning and trimming germline sequences using the consensus sequences. A cost of Inf will result in a gapless alignment.
  
• Default: 10
  
  

gap.extension.cost

• Description: The cost for extending a gap in when aligning and trimming germline sequences. A cost of Inf will result in a gapless alignment.
  
• Default: 4
  
  

The Af_plot_tree() function retrieves the igraph object from the provided AntibodyForests object for the specified clone within the specified sample and plots the lineage tree using the specified plotting parameters. This visualization offers insights into sequence evolution, enabling the observation of evolutionary processes in relation to sequence distances, clonal expansion, and other attributes such as isotype or evolutionary likelihood, as well as additional features appended to the VDJ dataframe.

To investigate the influence of different tree constructions methods on the resulting tree topologies, the AntibodyForests package allows for the comparison of trees generated using different algorithms. The Af_compare_methods() function can be used to compare the tree topologies of multiple AntibodyForests objects. This function calculates the distance between the trees of the different objects, clusters the trees, and visualizes the results. Various distance metrics can be used, such as the euclidean distance between the germline-to-node distance for each node or the Generalized Branch Length Distance (GBLD) between the trees. Visualization of the results can be done using a heatmap, an MDS or PCA plot. The function returns a list of the calculated distances, the clustering, and the visualization for each clonotype. Additionally the results can be averaged over all clonotypes.

When multiple AntibodyForests objects are created from the same VDJ dataframe and using the same sequence columns, the nodes in the resulting igraph objects should be the same. However, selecting different sequence columns may result in a different number of nodes and different node names. Additionally, the IgPhyML tool does not label the nodes according to their size but rather by using barcodes. When comparing tree topologies, having synchronized node names is essential.To achieve this, the Af_sync_nodes() function can be used. This function takes one AntibodyForests object as a reference and renames the nodes of all clonotypes within all samples of another AntibodyForests object to match the reference. The function matches the nodes using barcodes. If the barcodes of one node in the reference object are found in multiple nodes in the subject object, the nodes in the subject will get suffixes A, B, etc. (e.g., node4 will become node4A and node4B). Conversely, if the barcodes of multiple nodes in the reference object are found in only one node in the subject object, the node in the subject object will get the concatenated names of the nodes in the reference object (e.g., node2 and node4 will become node2+4).

The complex process of somatic hypermutation in a repertoire can be analyzed by comparing lineage tree topologies between clonotypes to identify different evolutionary patterns. The AntibodyForests package provides a function, Af_compare_within_repertoires(), to calculates the distance between all trees in the repertoire. This can be the euclidean distance between various calculated metrics, or the Jensen-Shannon distance between spectral density profiles of trees. The function can cluster the trees and visualize the results using a heatmap, MDS, or PCA plot. The function returns a list containing the distance matrix, the assigned clusters, and the visualization. The functions Af_cluster_metrics() and Af_cluster_node_features() can be used to analyze the differences between the clusters of trees.

The function Af_metrics() is used internally in different AntibodyForests function, but can also be used as standalone function to calculate metrics for each clonotype. Basic metrics such as the average germline-to-node depth, total number of nodes and cells, and sackin index can be calculated, as well as more complex metrics based on the Laplacian Spectrum of the trees. These metrics can be returned in a dataframe, or added to the AntibodyForests object.

The functions Af_distance_boxplot() and Af_distance_scatterplot can be used to observe patterns between the distance to the germline and various node features in the AntibodyForests object.

The quantify different evolutionary patterns between (groups of) repertoires within or across individuals, the AntibodyForests package provides a function, Af_compare_across_repertoires(), to calculates tree topology metrics of trees and compares them between repertoires. The function included different graph theory metrics such as the edge betweenness or the degree centrality, but it also includes all the metrics available in the Af_metrics() function. The function can plot these metrics in boxplots or frequency polygons per repertoire and calculates the statistical significance.

PLMs have demonstrated success in understanding the unique features of antibody evolution, such as germline gene usage and somatic hypermutation (SHM). The probabilities resulting from PLMs can be used as an estimate for evolutionary likelihood, both at the per-residue as the per-sequence (pseudolikelihood) level. These likelihoods can be incorporated into AntibodyForests to evaluate patterns in pseudolikelihood and tree topology, and to evaluate mutations along the edges of the trees.

The pseudolikelihoods and probability matrices can be supplied manually, or easily generated using the PLM-pipeline. Here, a VDJ-dataframe can be used as input CSV file, and the pipeline can generate the pseudolikelihoods and probability matrices of those sequences using different PLMs. The resulting pseudolikelihoods can be incoporates as node feature in the AntibodyForests-object using the Af_add_node_feature() function. The pseudolikelihoods can be visualized using the Af_plot_tree() function, and the correlation with distance to the germline can be visualized using the Af_distance_scatterplot function.

Sequences corresponding to nodes in the AntibodyForests-object can be extracted using the Af_get_sequences() function. This dataframe can be used as input CSV file for the PLM-pipeline to be able to match probability matrices to the right trees and nodes. PLM-generated probability matrices can be supplied to the function Af_PLM_dataframe() to generate a dataframe with the ranks and probabilities of the mutations along the edges of the lineage trees. These ranks and probabilties of both the mutated amino acids and the original amino acid can be visualized using the Af_plot_PLM() function.

AntibodyForests lineage trees reflect sequence-based variation during affinity maturation, but the 3D structure of the antibodies and antigens determine the actual binding. The 3D structure of antibodies can be predicted using computational methods such as homology modelling, or the structure can be inferred experimentally. AntibodyForests provides functions to incorporate 3D coordinates into the lineage trees, which allows the analysis of the evolution of the 3D structure together with somatic hypermutation.

AntibodyForests needs as input superimposed PDB files of protein structures corresponding to the sequences of the nodes in the trees. For example, the structures can be inferred with AlphaFold3 and superimposed on the C-alpha backbone using the Bio.PDB package in Python. Next to the PDB file, AntibodyForests can incorporate the 3Di sequence of the structures, which can be obtained using mini3di. Additionally, pKa values and free energy of antibody-antigen complexes can be integrated in the AntibodyForests object using the output of propka.

The function VDJ_3d_properties() can be used to calculate various biophysical properties from the PDB files, such as the charge, hydrophobicity, pLDDT, RMSD or 3Di-distance to the germline, or the pKa shift and free energy. These properties can be added to the AntibodyForests object using the Af_add_node_feature() function. The 3D properties can be visualized using the Af_plot_tree() function, and the correlation with distance to the germline can be visualized using the Af_distance_scatterplot function. The function Af_edge_RMSD() can be used to calculate the RMSD between the 3D structures of the mutating sequences along the edges of the lineage trees.